C  MODPATH release: Version 4.00 (V4, Release 1, 2-2000)
C    Changes to work with MODFLOW-2000
C
C MODPATH Version 3.00 (V3, Release 2, 5-99)
C Changes:
C   No change from previous release: (V3, Release 1, 9-94)
C***** SUBROUTINES *****
C     DRIVER
C***********************
 
C***** SUBROUTINE *****
      SUBROUTINE DRIVER (QX,QY,QZ,POR,IBOUND,XMAX,DELR,YMAX,DELC,ZBOT,
     2ZTOP,HEAD,BUFF,LAYCON,NCON,IBUFF,IUNIT,XLC,YLC,ZLC,
     3ZLLC,TOT,JLC,ILC,KLC,QSS,QSTO,
     4INILOC,PERLEN,NUMTS,TIMX,IBSTRT,HDRY,
     5HNOFLO,VER,LAYCBD)
C
      INCLUDE 'idat1.inc'
      COMMON /IDAT2/ IRCHTP,IEVTTP,NRCHOP,NEVTOP

      DIMENSION LAYCBD(200)
      DIMENSION TOT(NPART),XLC(NPART),YLC(NPART),ZLC(NPART),
     1ZLLC(NPART),JLC(NPART),ILC(NPART),KLC(NPART),QX(NCP1,NROW,NLAY),
     2QY(NCOL,NRP1,NLAY),QZ(NCOL,NROW,NLP1),POR(NCOL,NROW,NLPOR),
     3IBOUND(NCOL,NROW,NLAY),XMAX(NCOL),DELR(NCOL),YMAX(NROW),
     4DELC(NROW),ZBOT(NZDIM),ZTOP(NZDIM),
     5IBUFF(NCOL,NROW,NLAY),IUNIT(NUNIT),HEAD(NCOL,NROW,NLAY),
     6LAYCON(NLAY),NCON(NLAY),BUFF(NCOL,NROW,NLAY),QSS(NCOL,NROW,NLAY),
     7IBSTRT(NCOL,NROW,NLAY),
     9QSTO(NCOL,NROW,NLAY),INILOC(NPART)
      DIMENSION PERLEN(NPER),NUMTS(NPER),TIMX(NPER)
      CHARACTER*78 MES
      CHARACTER*(*) VER
C
C  GET TRIMMED LENGTH OF VERSION HEADER STRING
      CALL CHOP(VER,LVER)
 
C
      CALL DATIN (LAYCON,NCON,HEAD,DELT,XMAX,YMAX,DELR,DELC,
     1ZTOP,ZBOT,IBOUND,POR,QX,QY,QZ,IUNIT,IRP,NPRT,FRAC,BUFF,
     2IREV,IBUFF,JLC,ILC,KLC,XLC,YLC,ZLLC,ITMS,QSS,QSTO,ISNK,ITREAD,
     3IZSTOP,IPEP,MAXSTP,
     4INILOC,TSTOP,TIMREL,KKPER,KKSTP,PERLEN,NUMTS,
     5TIMX,IBSTRT,ZLC,HDRY,HNOFLO,TOT,ICMPCT,TBGABS,LAYCBD)
C
C  CHECK MASS BALANCES
C
      MES= 'DO YOU WANT TO COMPUTE VOLUMETRIC BUDGETS FOR ALL CELLS ?'
      CALL YESNO (MES,IANS,'3.1.4')
      IF(IANS.EQ.1) THEN
      CALL MPRMPT('SPECIFY AN ERROR TOLERANCE (IN PERCENT):')
      CALL GETR(1,ETOL,1,0.,0.,'3.1.1')
      IF(ISS.EQ.0) THEN
      NSTEP=0
      CALL READHQ(QX,QY,QZ,QSS,QSTO,HEAD,IBOUND,BUFF,NCOL,
     1 NROW,NLAY,NCP1,NRP1,NLP1,I10,NSTEP,DTSTP,TOTSIM,IERR,
     2 NUMTS,NPER,NSTEPF,NSTEPL,KPER,KSTP,0,I7)
      WRITE(I7,*) ' '
      WRITE(I7,9300) NSTEPF,NSTEPL
9300  FORMAT(1X,'COMPUTE BUDGETS FOR MODFLOW CUMULATIVE TIME STEPS ',I6,
     1' THROUGH ',I6)
      DO 5 NSTEP=NSTEPF,NSTEPL
      CALL READHQ(QX,QY,QZ,QSS,QSTO,HEAD,IBOUND,BUFF,NCOL,
     1 NROW,NLAY,NCP1,NRP1,NLP1,I10,NSTEP,DTSTP,TOTSIM,IERR,
     2 NUMTS,NPER,NSTEPF,NSTEPL,KPER,KSTP,1,I7)
      CALL BALNCE (QX,QY,QZ,QSS,QSTO,IBOUND,HEAD,NCOL,NROW,NLAY,
     1NCP1,NRP1,NLP1,NUMERR,I7,ISILNT,ETOL,KPER,KSTP,ISS,HDRY)
5     CONTINUE
      ELSE IF(ISS.EQ.1) THEN
      CALL BALNCE (QX,QY,QZ,QSS,QSTO,IBOUND,HEAD,NCOL,NROW,NLAY,
     1NCP1,NRP1,NLP1,NUMERR,I7,ISILNT,ETOL,1,1,ISS,HDRY)
      END IF
      END IF
C
C  CHECK DATA CELL BY CELL
C
      MES= ' DO YOU WANT TO CHECK DATA CELL BY CELL ?'
      CALL YESNO (MES,IANS,'3.1.2')
      IF(IANS.EQ.1) THEN
      CALL CELDAT (XMAX,YMAX,DELR,DELC,ZBOT,ZTOP,DELZCB,POR,IBOUND,HEAD,
     1 QX,QY,QZ,LAYCON,NCOL,NROW,NLAY,NZDIM,IGRID,NLPOR,NCP1,NRP1,NLP1,
     2 NCON,QSS,IBATCH,I7,QSTO,ISS,BUFF,NUMTS,NPER,I10,IBSTRT)
      END IF
C
      MES= 'SUMMARIZE FINAL STATUS OF PARTICLES IN SUMMARY.PTH FILE ?'
      CALL YESNO(MES,IOP,'3.1.3')
C
      IF (IBATCH.EQ.0) THEN
      WRITE (*,*) ' '
      WRITE (*,*) 'NOW COMPUTING PARTICLE PATHS...'
      WRITE (*,*) ' '
      IF(NPRT.EQ.0) THEN
        WRITE(*,*) 'NUMBER OF PARTICLES = 0. RUN STOPPED.'
		call stopfile  ! emrl
        STOP
      END IF
      END IF
 
      IF(NPRT.EQ.0) THEN
        WRITE(I7,*) 'NUMBER OF PARTICLES = 0. RUN STOPPED.'
		call stopfile  ! emrl
        STOP
      END IF
C
C
C... initialize for first pass through time loop
C
C... ISTOP is a flag that is used to indicate if another pass through
C       the time loop should be made.
C       ISTOP =  0 --> make another pass through time loop
C       ISTOP = -1 --> time at end of current modflow time step is
C                      greater than a user-specified stopping time
C                      (a negative value of ISTOP may be reset to 0
C                       if the current request time is less than the
C                       user-specified stop time)
C       ISTOP =  1 --> do not make another pass through the time loop
C
      ISTOP=0
 
      KOUNT=1
      TOLD=0.0E+0
      TIMREQ=0.0
 
      IF(ITREAD.EQ.1) READ(I8,*) MAXSTP,TFAC
      CALL GETREQ(ITMS,ITREAD,MAXSTP,KOUNT,DELT,TFAC,TIMREQ,I8)
      IF(ISS.EQ.0) THEN
      CALL GETSTP(KKPER,KKSTP,NUMTS,NPER,I7,NSTEP)
      CALL READHQ(QX,QY,QZ,QSS,QSTO,HEAD,IBOUND,BUFF,NCOL,
     1 NROW,NLAY,NCP1,NRP1,NLP1,I10,NSTEP,DTSTP,TOTSIM,IERR,
     2 NUMTS,NPER,NSTEPF,NSTEPL,KPER,KSTP,1,I7)
      IF(IERR.EQ.1) THEN
      WRITE(*,*) 'TIME STEP NOT FOUND IN COMPOSITE BUDGET FILE. STOP.'
      WRITE(I7,*) 'TIME STEP NOT FOUND IN COMPOSITE BUDGET FILE. STOP.'
		call stopfile  ! emrl
      STOP
      ELSE
      WRITE(I7,*) ' '
      WRITE(I7,9000) KKPER,KKSTP,NSTEP
9000  FORMAT(1X,'INITIAL BUDGET AND HEAD DATA READ FROM CBF FOR:'
     1/10X,'PERIOD ',I4,2X,'STEP ',I5,2X,'(CUMULATIVE TIME STEP = ',I6,
     2')')
      END IF
      IF(IREV.EQ.1) THEN
      TIMSTP=TIMREL*DTSTP
      ELSE
      TIMSTP= (1.0-TIMREL)*DTSTP
      END IF
      ELSE IF(ISS.EQ.1) THEN
      NSTEP=1
      TIMSTP=1.0E+30
      END IF
      IF(IREV.EQ.1) CALL FLIPQ(QX,QY,QZ,NCOL,NROW,NLAY,NCP1,NRP1,NLP1)
 
      IF(TSTOP.LE.TIMSTP) THEN
      ISTOP= -1
      TIMSTP=TSTOP
      END IF
 
C  WRITE VERSION HEADER STRING ON LINE 1 OF TIMESERIES FILE
      IF(MODE.EQ.2) CALL FILHDR(I4,VER,ICMPCT,TBGABS)
 
C
C  COMPUTE GLOBAL Z COORDINATES OF PARTICLES THAT ARE ACTIVE AT TIME 0.
C  IF TIME SERIES MODE, WRITE INITIAL CONDITIONS IN TIMESERIES FILE
C
        DO 50 N=1,NPRT
          IF(TOT(N).EQ.-1.0) THEN
            CALL ZCALC(IGRID,NCOL,NROW,NLAY,ZLC(N),ZLLC(N),ZTOP,ZBOT,
     1               HEAD(JLC(N),ILC(N),KLC(N)),LAYCON(KLC(N)),
     2               ILC(N),JLC(N),KLC(N),NZDIM,I7)
            READ(I5,REC=N) XSTRT,YSTRT,ZLSTRT,IPCODE,TRLEAS,NDUMMY
            WRITE(I5,REC=N) XSTRT,YSTRT,ZLSTRT,IPCODE,TRLEAS,NSTEP
            IF (MODE.EQ.2) THEN
              CALL WRITTS(I4,ICMPCT,0,N,JLC(N),ILC(N),KLC(N),XLC(N),
     1          YLC(N),ZLC(N),ZLLC(N),0.0,NSTEP,NROW,NCOL)
            END IF
          END IF
50      CONTINUE
 
C
C  FIND THE TIME OUT CONDITION FOR THE FIRST PASS THROUGH THE TIME LOOP
C
C... ICASE = 1 --> request time > end of current modflow time step
C                  (new flows must be loaded at the end of this pass
C                   through the time loop)
C... ICASE = 2 --> request time < end of current modflow time step
C                  (existing flows can be used for this pass and the
C                   next pass through the time loop)
C... ICASE = 3 --> request time = end of current modflow time step
C                  (new flows must be loaded at the end of this pass
C                   through the time loop)
      IF(TIMREQ.GT.TIMSTP) THEN
      ICASE=1
      TMAX=TIMSTP
      ELSE IF(TIMREQ.LT.TIMSTP) THEN
      ICASE=2
      TMAX=TIMREQ
      IF(ISTOP.LT.0) ISTOP=0
      ELSE IF(TIMREQ.EQ.TIMSTP) THEN
      ICASE=3
      TMAX=TIMSTP
      END IF
 
C... WRITE VERSION HEADER STRING TO PATHLINE FILE
      IF(MODE.EQ.1) CALL FILHDR(I2,VER,ICMPCT,TBGABS)
 
C
C---------------------  START TIME LOOP  -------------------------
C
C... imprint the zone codes onto the IBOUND array
      DO 61 K=1,NLAY
      DO 61 I=1,NROW
      DO 61 J=1,NCOL
61    IBOUND(J,I,K)= IBSTRT(J,I,K)*IBOUND(J,I,K)
 
 
60    CONTINUE
 
C
C--------------------  START PARTICLE LOOP  -----------------------
C
      IEND=0
      DO 70 N=1,NPRT
C
C  SET START TIME FOR PARTICLE
C
      TIME=TOLD
 
C  IF RUN IS FORWARD MODE, CHECK TO SEE IF PARTICLE IS ACTIVE OR READY
C  TO BE RELEASED.
      IF(IREV.EQ.0) THEN
        IF(TOT(N).EQ.-1.0E+30) THEN
          NODE=INILOC(N)
          CALL RELEAS(N,TOT(N),ZTOP,ZBOT,HEAD,LAYCON,ZLC(N),ZLLC(N),
     1         NODE,NZDIM,NCOL,NROW,NLAY,TIME,TMAX,I5,I7,IGRID,
     2         HDRY,HNOFLO,IPSTAT,NSTEP)
        END IF
      END IF
 
C  IF PARTICLE HAS DISCHARGED, SKIP TO THE NEXT PARTICLE
      IF(TOT(N).GE.0.0E+0) GO TO 70
 
C  IF PARTICLE IS UNRELEASED, SKIP TO THE NEXT PARTICLE
      IF(TOT(N).EQ.-1.0E+30) THEN
        IF(IPSTAT.EQ.1) IEND=IEND+1
        GO TO 70
      END IF
 
C  IF IT GETS THIS FAR, THEN THERE IS AT LEAST 1 ACTIVE PARTICLE TO TRACK
C  THROUGH THIS PASS THROUGH THE TIME LOOP. SET NSTEPC=NSTEP TO MAKE SURE
C  THAT NSTEPC WILL CONTAIN THE VALUE THAT NSTEP HAD DURING THE FINAL PASS
C  THROUGH THE TIME LOOP
C
      NSTEPC=NSTEP
 
C  MOVE PARTICLE AS FAR AS IT CAN GO THIS PASS THROUGH THE TIME LOOP
C
      CALL FLOLIN (N,IGRID,MODE,TIME,TMAX,IDSCH,
     1JLC(N),ILC(N),KLC(N),XLC(N),YLC(N),ZLC(N),ZLLC(N),IBOUND,LAYCON,
     2ZBOT,ZTOP,XMAX,YMAX,QX,QY,QZ,DELC,DELR,POR,HEAD,NCON,NCOL,NROW,
     3NLAY,NLPOR,NZDIM,NCP1,NRP1,NLP1,ISNK,IREV,FRAC,IZSTOP,I2,
     4I7,QSS,QSTO,INILOC,NPART,ISS,ICASE,HDRY,NSTEP,ICMPCT)
 
C  ASSIGN TIME OF TRAVEL TO TOT ARRAY IF PARTICLE HAS DISCHARGED
C
C  IDSCH = -2 -- PARTICLE WAS STOPPED IN AUTOMATIC DISCHARGE ZONE
C  IDSCH = -1 -- PARTICLE WAS STRANDED IN INACTIVE CELL
C  IDSCH =  0 -- PARTICLE HAS DISCHARGED NORMALLY
C  IDSCH =  1 -- PARTICLE HAS NOT DISCHARGED
C
      IF(IDSCH.LE.0) THEN
      TOT(N)=TIME
      READ(I5,REC=N) XSTRT,YSTRT,ZLSTRT,IPCODE,TRLEAS,NDUMMY
      IF(IDSCH.EQ.0) IPCODE=  1 + (10*NSTEP)
      IF(IDSCH.EQ.-1) IPCODE= -1
      IF(IDSCH.EQ.-2) IPCODE= 2 + (10*NSTEP)
      WRITE(I5,REC=N) XSTRT,YSTRT,ZLSTRT,IPCODE,TRLEAS,NDUMMY
      ELSE
      IEND=IEND+1
      END IF
C
C  DON'T WRITE RESULTS IN TIMESERIES FILE FOR TIME=TMAX IF PARTICLE
C  DISCHARGED AT TIME < TMAX.
C
      IF(IDSCH.LE.0.AND.TIME.LT.TMAX) GO TO 70
C
C  WRITE RESULTS TO TIMESERIES FILE
C
      IF (MODE.EQ.2 .AND. ICASE.GT.1) THEN
        CALL WRITTS(I4,ICMPCT,KOUNT,N,JLC(N),ILC(N),KLC(N),XLC(N),
     1   YLC(N),ZLC(N),ZLLC(N),TMAX,NSTEP,NROW,NCOL)
      END IF
C
70    CONTINUE
C
      IF(IEND.EQ.0) GO TO 71
        IF(ISTOP.LT.0) THEN
          ISTOP=1
          GO TO 71
        END IF
 
C... if we are entering a new modflow time step on the next pass through the
C    time loop, then get new flow rates.
      IF(ICASE.EQ.1 .OR. ICASE.EQ.3) THEN
 
        IF(ISS.EQ.1) THEN
          ISTOP=1
          GO TO 71
        END IF
 
        IF(IREV.EQ.1) THEN
          NSTEP=NSTEP-1
        ELSE
          NSTEP=NSTEP+1
        END IF
 
C... STOP IF WE HAVE REACHED THE END OF THE CBF
        IF(NSTEP.LT.NSTEPF .OR. NSTEP.GT.NSTEPL) THEN
          ISTOP=1
          GO TO 71
        END IF
 
        CALL READHQ(QX,QY,QZ,QSS,QSTO,HEAD,IBOUND,BUFF,NCOL,
     1  NROW,NLAY,NCP1,NRP1,NLP1,I10,NSTEP,DTSTP,TOTSIM,IERR,
     2  NUMTS,NPER,NSTEPF,NSTEPL,KPER,KSTP,1,I7)
 
        IF(IERR.NE.0) THEN
          ISTOP=1
          WRITE(I7,9100) NSTEP
9100      FORMAT(1X,'ERROR READING COMPOSITE BUDGET FILE:',
     1           /10X,'CUMULATIVE TIME STEP = ',I6)
          GO TO 71
        END IF
        WRITE(I7,9200) NSTEP
9200    FORMAT(1X,'NEW MODFLOW TIME STEP: DATA READ FROM CBF; CUMULATIVE
     1 TIME STEP = ',I6)
 
        IF(IREV.EQ.1) THEN
          CALL FLIPQ(QX,QY,QZ,NCOL,NROW,NLAY,NCP1,NRP1,NLP1)
        END IF
 
        TIMSTP=TIMSTP+DTSTP
        IF(TSTOP.LE.TIMSTP) THEN
          ISTOP= -1
          TIMSTP=TSTOP
        END IF
 
C... IMPRINT THE ZONE CODES ONTO THE IBOUND ARRAY
      DO 62 K=1,NLAY
      DO 62 I=1,NROW
      DO 62 J=1,NCOL
62    IBOUND(J,I,K)= IBSTRT(J,I,K)*IBOUND(J,I,K)
 
      END IF
 
C... set the end time for the next pass through the time loop and
C    determine if it will be necessary to load new flows at the end
C    of the pass. get the new request time if necessary.
      IF(ICASE.EQ.2 .OR. ICASE.EQ.3) THEN
        KOUNT=KOUNT+1
        CALL GETREQ(ITMS,ITREAD,MAXSTP,KOUNT,DELT,TFAC,TIMREQ,I8)
      END IF
      TOLD=TMAX
      IF(TIMREQ.GT.TIMSTP) THEN
      ICASE=1
      TMAX=TIMSTP
      ELSE IF(TIMREQ.LT.TIMSTP) THEN
      ICASE=2
      TMAX=TIMREQ
      IF(ISTOP.LT.0) ISTOP=0
      ELSE IF(TIMREQ.EQ.TIMSTP) THEN
      ICASE=3
      TMAX=TIMSTP
      END IF
C
71    CONTINUE
C-------------------  END PARTICLE LOOP  -------------------------
C
C  CHECK TO SEE IF ALL PARTICLES HAVE DISCHARGED. IF NOT GO THROUGH
C  TIME LOOP AGAIN.
C
      IF(IEND.GT.0.AND.ISTOP.NE.1) GO TO 60
      IF(ISTOP.EQ.1) THEN
      DO 75 N=1,NPRT
      IF(TOT(N).LT.0.0 .AND. TOT(N).NE.-1.0E+30) TOT(N)= TMAX
75    CONTINUE
      END IF
C
C----------------------  END TIME LOOP  ---------------------------
C
C  WRITE DISCHARGE DATA TO ENDPOINT  FILE
C
      IACTIV=0
      ISTRAN=0
      ISTOPD=0
      IDCH=0
      IUNREL=0
 
C  WRITE VERSION HEADER AS FIRST LINE OF ENDPOINT FILE
      CALL FILHDR(I3,VER,ICMPCT,TBGABS)
 
      TIMMIN=1.0E+30
      TIMMAX=0.0
      KNTPTS=0
      TIMSUM=0.0
      DO 80 N=1,NPRT
      READ(I5,REC=N) XSTRT,YSTRT,ZLSTRT,IPCODE,TRLEAS,NSFRST
      IF(IPCODE.GE.0) THEN
        IDCODE= MOD(IPCODE,10)
      ELSE
        IDCODE= IPCODE
      END IF
      NODE=INILOC(N)
      IF(NODE.LT.0) NODE= -NODE
      CALL ND2IJK(NODE,IFRST,JFRST,KFRST,NROW,NCOL)
      IF(IDCODE.EQ.0) THEN
        IACTIV=IACTIV+1
      ELSE IF(IDCODE.EQ.1) THEN
        IDCH=IDCH+1
      ELSE IF(IDCODE.EQ.-1) THEN
        ISTRAN=ISTRAN+1
      ELSE IF(IDCODE.EQ.2) THEN
        ISTOPD=ISTOPD+1
      ELSE IF(IDCODE.EQ.-2) THEN
        IUNREL=IUNREL+1
      END IF
C
C  WRITE DISCHARGE INFO TO SUMMARY FILE FOR 1ST 1000 PARTICLES IF IOP=1
C
      IF(IOP.EQ.1) THEN
      IF(N.EQ.1) THEN
      WRITE(I7,*) ' '
      WRITE(I7,*)' ******** PARTICLE TERMINATION INFORMATION'
      END IF
      IF(N.LE.1000) THEN
        IF(IDCODE.EQ.2) THEN
          WRITE (I7,5010) N,ILC(N),JLC(N),KLC(N),TOT(N)
        ELSE IF(IDCODE.EQ.1) THEN
          WRITE (I7,5013) N,ILC(N),JLC(N),KLC(N),TOT(N)
        ELSE IF(IDCODE.EQ.-1) THEN
          WRITE (I7,5012) N,ILC(N),JLC(N),KLC(N),TOT(N)
        ELSE IF(IDCODE.EQ.0) THEN
          WRITE (I7,5011) N,ILC(N),JLC(N),KLC(N),TOT(N)
        ELSE IF(IDCODE.EQ.-2) THEN
          WRITE (I7,5014) N,ILC(N),JLC(N),KLC(N),TOT(N)
        END IF
      ELSE IF(N.EQ.1001) THEN
        WRITE(I7,*)
     1'***** DATA SUMMARIZED FOR FIRST 1000 PARTICLES ONLY ****'
        END IF
      END IF
C
      IZONE=IBOUND(JLC(N),ILC(N),KLC(N))
      IF(IZONE.LT.0) IZONE=-IZONE
      IF(IZONE.GE.1000) IZONE=IZONE/1000
      IZONE2=IBOUND(JFRST,IFRST,KFRST)
      IF(IZONE2.LT.0) IZONE2=-IZONE2
      IF(IZONE2.GE.1000) IZONE2=IZONE2/1000
 
C  WRITE ENDPOINT DATA FOR A PARTICLE
      IF(IPCODE.EQ.0) IPCODE= 10*NSTEPC
      IF (IPEP.EQ.0) THEN
        IF(TOT(N).GE.0.0 .AND. IPCODE.NE.-2) THEN
          KNTPTS=KNTPTS+1
          TIMSUM=TIMSUM+(TOT(N)-TRLEAS)
          IF(TOT(N)-TRLEAS.LT.TIMMIN) TIMMIN=TOT(N)-TRLEAS
          IF(TOT(N)-TRLEAS.GT.TIMMAX) TIMMAX=TOT(N)-TRLEAS
        END IF
        ISKIP=0
        IF(TOT(N).LT.0.0 .OR. IPCODE.EQ.-2) ISKIP=1
        IF(MODE.GT.0) ISKIP=0
        IF(ISKIP.EQ.0) THEN
        CALL WRITEP(I3,ICMPCT,IZONE,JLC(N),ILC(N),KLC(N),XLC(N),
     1   YLC(N),ZLLC(N),ZLC(N),TOT(N),XSTRT,YSTRT,ZLSTRT,JFRST,IFRST,
     2   KFRST,IZONE2,NSFRST,IPCODE,TRLEAS,NROW,NCOL)
        END IF
      ELSE IF(IZONE.EQ.IZSTOP) THEN
        IF(TOT(N).GE.0.0 .AND. IPCODE.NE.-2) THEN
          KNTPTS=KNTPTS+1
          TIMSUM=TIMSUM+(TOT(N)-TRLEAS)
          IF(TOT(N)-TRLEAS.LT.TIMMIN) TIMMIN=TOT(N)-TRLEAS
          IF(TOT(N)-TRLEAS.GT.TIMMAX) TIMMAX=TOT(N)-TRLEAS
        END IF
        ISKIP=0
        IF(TOT(N).LT.0.0 .OR. IPCODE.EQ.-2) ISKIP=1
        IF(MODE.GT.0) ISKIP=0
        IF(ISKIP.EQ.0) THEN
        CALL WRITEP(I3,ICMPCT,IZONE,JLC(N),ILC(N),KLC(N),XLC(N),
     1   YLC(N),ZLLC(N),ZLC(N),TOT(N),XSTRT,YSTRT,ZLSTRT,JFRST,IFRST,
     2   KFRST,IZONE2,NSFRST,IPCODE,TRLEAS,NROW,NCOL)
        END IF
      END IF
 
80    CONTINUE
      TIMAVE=0.
      IF(KNTPTS.GT.0) TIMAVE= TIMSUM/FLOAT(KNTPTS)
      KNTPT2=0
      DO 81 N=1,NPRT
        READ(I5,REC=N) XSTRT,YSTRT,ZLSTRT,IPCODE,TRLEAS,NSFRST
        IF(IPCODE.GE.0) THEN
          IDCODE= MOD(IPCODE,10)
        ELSE
          IDCODE= IPCODE
        END IF
        IF(IPEP.GT.0 .AND. IDCODE.NE.2) GO TO 81
        IF(TOT(N).GE.0.0 .AND. IPCODE.NE.-2) THEN
          IF(TOT(N)-TRLEAS.LE.TIMAVE) KNTPT2=KNTPT2+1
        END IF
81    CONTINUE
      TIMFR=0.0
      IF(KNTPT2.GT.0) TIMFR= 100.*FLOAT(KNTPT2)/FLOAT(KNTPTS)
      ISUM=IACTIV+IDCH+ISTOPD+ISTRAN
      IF(IREV.EQ.0) ISUM=ISUM+IUNREL
      IF(ISILNT.EQ.0) THEN
        WRITE(*,*)
        IF(IPEP.EQ.0) THEN
          WRITE(*,5036)
        ELSE IF(IPEP.GT.0) THEN
          WRITE(*,5037) IZSTOP
        END IF
        WRITE(*,5034) TIMMIN,TIMMAX,TIMAVE
        WRITE(*,5035) TIMFR
        WRITE(*,*)
        IF(IREV.EQ.0) THEN
          WRITE(*,5032) IACTIV,IDCH,ISTOPD,ISTRAN,IUNREL,ISUM,NPRT
        ELSE
          WRITE(*,5031) IACTIV,IDCH,ISTOPD,ISTRAN,ISUM,NPRT
        END IF
      END IF
 
      WRITE(I7,*) '   '
      IF(IPEP.EQ.0) THEN
        WRITE(I7,5036)
      ELSE IF(IPEP.GT.0) THEN
        WRITE(I7,5037) IZSTOP
      END IF
      WRITE(I7,5034) TIMMIN,TIMMAX,TIMAVE
      WRITE(I7,5035) TIMFR
      WRITE(I7,*) '   '
      IF(IREV.EQ.0) THEN
        WRITE(I7,5032) IACTIV,IDCH,ISTOPD,ISTRAN,IUNREL,ISUM,NPRT
      ELSE
        WRITE(I7,5031) IACTIV,IDCH,ISTOPD,ISTRAN,ISUM,NPRT
      END IF
 
5010  FORMAT(' PARTICLE ',I6,'  WAS STOPPED IN ROW ',I3,', COL ',I3,
     1', LAY ',I3,'  TIME= ',1PE11.4)
5011  FORMAT(' PARTICLE ',I6,'    IS ACTIVE IN ROW ',I3,', COL ',I3,
     1', LAY ',I3,'  TIME= ',1PE11.4)
5012  FORMAT(' PARTICLE ',I6,' WAS STRANDED IN ROW ',I3,', COL ',I3,
     1', LAY ',I3,'  TIME= ',1PE11.4)
5013  FORMAT(' PARTICLE ',I6,'        EXITS IN ROW ',I3,', COL ',I3,
     1', LAY ',I3,'  TIME= ',1PE11.4)
5014  FORMAT(' PARTICLE ',I6,' WAS NOT RELEASED')
5030  FORMAT(4(1X,I4),8(1X,E13.6),4(1X,I4),I3,1X,E13.6)
5031  FORMAT(1X,I7,' PARTICLES REMAIN ACTIVE'/
     11X,I7,' PARTICLES STOPPED AT INTERNAL SINKS/SOURCES OR BOUNDARIES'
     2/1X,I7,' PARTICLES STOPPED IN AN AUTOMATIC TERMINATION ZONE'/
     31X,I7,' PARTICLES WERE STRANDED IN INACTIVE CELLS'/
     41X,'-----'/1X,I7,' PARTICLES ACCOUNTED FOR OUT OF A TOTAL OF ',I7)
5032  FORMAT(1X,I7,' PARTICLES REMAIN ACTIVE'/
     11X,I7,' PARTICLES STOPPED AT INTERNAL SINKS/SOURCES OR BOUNDARIES'
     2/1X,I7,' PARTICLES STOPPED IN AN AUTOMATIC TERMINATION ZONE'/
     31X,I7,' PARTICLES WERE STRANDED IN INACTIVE CELLS'/
     41X,I7,' PARTICLES WERE NOT RELEASED'/
     41X,'-----'/1X,I7,' PARTICLES ACCOUNTED FOR OUT OF A TOTAL OF ',I7)
5034  FORMAT(1X,'MINIMUM TRAVEL TIME = ',1PE13.5/1X,'MAXIMUM TRAVEL TIME
     1 = ',E13.5/1X,'AVERAGE TRAVEL TIME = ',E13.5)
5035  FORMAT(1X,F6.1,'% OF THE PARTICLES HAD TRAVEL TIMES LESS THAN THE
     1AVERAGE TRAVEL TIME')
5036  FORMAT(1X,'TRAVEL TIME SUMMARY FOR ALL PARTICLES:')
5037  FORMAT(1X,'TRAVEL TIME SUMMARY FOR PARTICLES TERMINATING IN ZONE',
     11X,I4,':')
C
6000  FORMAT(A4,A,A2)
      RETURN
      END
